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of distinct fourth order Diffusion Monte Carlo algorithms. These sophisticated 
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algorithms require higher derivatives of the drift velocity and local energy 
and are more complicated to program. However, they allowed very large time 
steps to be used, converged faster with lesser correlations, and virtually elim 
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lution operator to fourth order with positive coefficients, we derived a number 



inated the step size error. We demonstrated the effectiveness of these quartic 



algorithms by solving for the ground state properties of bulk liquid helium. 
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I. INTRODUCTION 

The basic idea of the Diffusion Monte Carlo (DMC) algorithm is to solve for the ground 
state of the Hamiltonian H by evolving the imaginary time Schrodinger equation 



-— V>(x,t) = i^(x,£) 



-V 2 + V(x) V(x,t) 



(1) 



to large time |T|-§[. Here, x and V 2 denote the coordinate and the Laplacian of the N- 
particle system. In order for the algorithm to be practical, capable of handling rapidly 
varying potentials, it is essential to implement important sampling as suggested by Kalos 
Q. This means that instead of solving for ^(x), one evolves the product wave function 
p(x) = 0(x) , 0(x) according to @,||] 



--p(x,t)=0(x)F</»- 1 (x)p(x,t), 
1, 



V 2 p(x, t) + Vi Gi(x)p(x, t) + £ L (x)p(x, t) 

z L J 

T + D + E L ]p(x,t) = Hp(x,t), 



(2) 
(3) 



where 



is the local energy 



£ L (x) = <j>(x)~ l H<p(x), 



(4) 



G<(x) = ^(x^Vi^x) = -V^(x) 



(5) 



is the <in/£ velocity and </>(x) = exp[— 5(x)] is the trial ground state wave function. 
Eq.(H) has the formal operator solution 



p{t) = e 



-t(T+D+E L ) 



P(0) 



-e(T+D+E L ) 



p(0). 



(6) 



Various DMC algorithms correspond to different approximations to the short time evolution 
operator e - £ ( T + D + B i) . Initial implementations [0-0] of the DMC algorithm correspond to 
essentially approximating 



e -e(T+D+E L ) ^ e -l € E Le -eT e -eD e -±eE L ^ Sj\ 

which is at most first order in e. By use of various clever tricks, this error can be reduced 
substantially in specific applications ||. However, it was recognized by Chin || that in 
order to have a general second order DMC algorithm, one must simulate the embedded 
Fokker-Planck evolution operator e~ e ^ T+D \ i.e. the Fokker-Planck equation 

- -p(x, t) = -^V 2 p(x, t) + Vi [Gf i (x)p(x, t)] = Lp(x, t), (8) 

correctly to second order. The reason for this is clear. In the limit when the trial function 
is the exact ground state wave function 0(x) — * ^o( x ), the local energy is the exact ground 
state energy, which is just a constant. The convergence of the DMC algorithm would then 
coincide with the convergence of the Langevin algorithm for simulating the Fokker-Planck 
equation. Thus in order to have a second order DMC algorithm |J, one must have a second 
order Langevin algorithm, for example, by approximating 

e -eL = e -e(T + D) ^ ^eT^eD^eT (g) 

This idea of operator factorization seemed promising for generating higher order DMC 
algorithms. However, Suzuki proved in 1991 that, beyond second order, it is impossible 
to factorize 

N 

exp[e(A + B)\ = JJ expfa^eA] exp^eS] (10) 

without having some coefficients Oj and hi being negative. Since e~ at€T is the diffusion 
kernel, a negative a% would imply a diffusion process backward in time, which is impossible 
to simulate. Thus higher than second order DMC algorithms cannot be based on obvious 
factorizations of the form (|10|). 

In this work, we show how to derive a number of distinct quartic DMC algorithms by 
factorizing the operator e - e ( T + D + E L) ^ fourth order with positive coefficients. We first 
review how each factorized operator can be simulated in Section II, followed by a derivation 
of a fourth order DMC algorithm in Section III. The backbone of this algorithm is a 4th 



order Langevin algorithm, which is of importance in its own right. In Section IV we examine 
the working details of this algorithm and check its quartic convergence on various systems 
including the practical case of liquid Helium. In Section V we discuss alternative quartic 
algorithms by considering the unrestricted factorization of e - e ( T + D + E L) _ ^^g convergences 
of two alternative 4th algorithms are also tested on liquid helium. Our conclusions and 
suggestions for future work are contained in Section VI. 

II. SIMULATING THE BASIC OPERATORS 

The method of operator factorization depends on the fact that each component factor 
can be simulated exactly or to the required order. The effect of e~ eT on p(x, t) is to evolve 
the latter forward in time according to the diffusion equation 

-|p(x,t) = -^V 2 p(x,t). (11) 

For a set of points {x{\ distributed according to p(x, £), this can be exactly simulated by 
updating each point according to 

x'i = Xi + v 7 ^, (12) 

where {£j} is a set of Gaussian distributed random numbers with zero mean and unit vari- 
ance. The operator e~ eD evolves p(x, t) forward in time according to the continuity equation 

-^p(x,t) = ^[G l (x)p(x,t)], (13) 

where Gj(x)p(x, t) = Jj(x) is the particle current density with drift velocity field Gj(x). 
This can also be exactly simulated by setting 

x[ = Xi(e), (14) 

where Xi(e) is the exact trajectory determined by 



d i - g «- < i5 > 



with the initial condition Xj(0) = Xj. In practice, one can only solve this trajectory equa- 
tion to the required order of accuracy. The operator e~ eEh evolves p(x, t) forward in time 
according to the rate equation 

--p(x,t)= J E; L (x)p(x,t). (16) 

The exact solution 

p(x,t + e) =e- e£; ^ x V(x,t) (17) 

can be simulated by updating the weight W k associated with the configuration x^ by 

W' k = e~ e[EL{xk) - E] W k . (18) 

A uniform constant E is usually added to keep the weights near unity. 

There are various methods |8]-[rTl| of keeping track of weights, the original and simplest 
method || is just to replicate the configuration x» on the average e - e i E L(^)~E] times. We 



use a method which is intermediate between that of || and [T(| • Our algorithm is, however, 
independent of any specific method of weight tracking. 
Thus the second order factorization of 

leads to the following second order DMC algorithm: 

Vi = Xi + Zi^Je/2, 

x >= yi (e)+C^/2, (20) 

where the final position x' is to be weighted by 

W = e -±e[E L (x)+E L (x')-2E]^ ^1) 

and where each yi(e) needs to be solved at least to second order by any convenient method. 
Following H this algorithm will be referred to as DMC2b. 



III. A FOURTH ORDER ALGORITHM 

For a fourth order factorization of exp[e(A + B)} with positive coefficients, Suzuki ]7| has 
shown that it is necessary to retain as a factor, the exponential of either double commutators 
[A, [B,A]] or [B, [A, B]]. Recently, Chin |E| has derived three such factorization schemes, 
two of which were also found previously by Suzuki [I3|. To decompose e - e ( T + D + E L) — 



e 6 ( l+e l) to fourth order, one possibility is to keep the Langevin operator L intact. In this 
case, the double commutator 

[E L , [L, E L \] = [E L , [T, E l \] = (diE L )(diE L ), (22) 

is the square of the gradient of the local energy, which is a manageable coordinate function. 
Since the Langevin operator is complicated to simulate, we must choose a fourth order 
factorization of e~ e< y L+EL ) which minimizes the appearance of L. We choose the following 



factorization as given by Refs. [12,13 



e -e(L+E L ) = e -leE Le -leL e -leE Le -±eL e — 6 eE L + Q ^ ^ 

with E L given by 

E L = E L + ^e 2 [E L , [L, E l }} = E L + ^e 2 1 VE L \ 2 . (24) 

Thus to the extend that the local energy -El(x) is a smooth function, the double commutator 
correction will be negligible. 

The weights in fl2lf ) have a simple structure. If x is the initial configuration, Xi/ 2 the 
Langevin evolved configuration time step e/2 later, and xi the Langevin evolved configura- 
tion a time step e/2 later still, then we assign the final configuration x x a weight of 

W ± = e -t[l E L(xi)+lE L (x 1/2 )+±E L (xo)-E}_ ( 25 ) 

The demanding part of this DMC algorithm is the simulation the Fokker-Planck equation 
(^). The resulting Langevin algorithm is an important simulation algorithm with numerous 
applications in statistical and chemical physics |14| . Since we have recently given a detailed 



6 



derivation of a fourth order Langevin algorithm [plj , we will be brief here in summarizing its 
essential features. To obtain a fourth order Langevin algorithm, we again seek to decompose 
e~ tL = e~ t( - T+D ^ to fourth order. In this case, we keep the double commutator [D, [T,D]], 
which is at most a second order differential operator, and factorize the Fokker-Planck oper- 
ator as JT2| 1, 

l^-^T^-hD-^f^UD -1(1 i) £ T 



e -eL = e -e(T+D) = e ~ 



V3^ e -2^ e V-i e -2^ e 



^ J "+o(e 5 ), 



(26) 



^f = -jf + E -{2 - V3)[D, [T, D]] = ^f+ € -{2 - VsM^f^ + 
where subscripts indicate partial differentiations, and 

J i>3 = ^i,k^j,k ^i,j,k^k 



OiVi 



(27) 



(28) 
(29) 



By appropriate normal ordering, the double commutator term can be regarded as a non- 
uniform Gaussian random walk. However, in order to be able to sample the non-uniform 
Gaussian in cases where fij has negative eigenvalues, we implement the normal ordering as 
follows so that the full covariance matrix is always positive definite in the limit of small e, 

^3 



exp W, 



cxp 



- (2 - Vs) (aa,/y + di 



AMexp 



tef)4 3 

-didj8i/j + |^ (2 - V3) {dtdjfij + diVi) \ exp f ^ T J > 



(30) 



_2^\ 2 

where jV denotes the normal ordering of all derivative operators to the left. Factorization 
(EST) can now be simulated as 



Wi = Xi + £ 
y l = w t (e/2)+£' l 



\^\ l A 



2^3' 



Zi =yi 



4=Zi{£/2)+ ek{i-±), 



^ + iS-r^» 



(,■ 



.1 • 



(31) 



where £» to £f are four sets of independent Gaussian random numbers with zero mean and 
unit variance. Here, the two trajectory equations Wi(e/2), Zi(e/2) must be solved correctly 
to at least fourth order. Empirically one observes that the more accurately one solves the 
trajectory equation, the smaller is the fourth order error coefficient. However, in practice 
one must weight improved convergence, which allows larger time steps to be used, against 
greater computational effort. In the present case, we solve the trajectory equation by the 
standard 4th order Runge-Kutta algorithm. 

Eq. fl23|) is our basic 4th order DMC algorithm and will be referred to as DM C4. We will 
first explore its workings in some detail before considering alternative algorithms. 

IV. APPLYING THE FOURTH ORDER ALGORITHM 

We begin by verifying that DMC4 is indeed quartic by solving the 3-D harmonic oscillator 

H = -iv 2 + ir 2 (32) 

both analytically with the help of Mathematica and by direct Monte Carlo simulation. The 
trial function used is 

0(r) = exp(--ar 2 ) (33) 

with a deliberate poor choice of the trial parameter a = 1.8. In Fig.l we plot the ground 
state energy from the mixed expectation 

E - ~m- (34) 

as a function of the step size e used. The lines are analytical functions from Mathematica and 
the plotting symbols are Monte Carlo simulation results. We have included one first order, 
two second order, and one first order rejection DMC algorithm for comparison. The detail 
descriptions of DMC1, DMC2a and DMC2b can be found in Ref. ||. The rejection algorithm 
uses a first order Langevin algorithm together with a generalized Metropolis acceptance- 
rejection step so that the square of the trial function is exactly sampled at all step sizes |§. 

8 



For this case, we have no analytical result and the plotted line is just a 6th order polynomial 
fit. The 4th order algorithm is indeed quartic, as can be verified analytically. It is distinct 
from lower order algorithms even in Monte Carlo simulations. 

We next test DMC4 by solving the 3-D Morse potential with Hamiltonian, 



H = --V 2 + D e [e" 2a ( r - r '°) - 2e- Q(r - ro) l , (35) 



with D e = 50, r = 1 and a = 10. These values ensure that the ground state is high up 
in the potential and that the ground state wave function is not well approximated by a 
Gaussian. We use a trial function of the form 

0(r) = exp(— ar — br~ 3 ), (36) 

with a = 15.29, b = 6.82 and variational energy -11.1774. This is to be compared with 
the exact ground state energy of -12.5 . The convergence of various DMC algorithms are 
compared in Fig.2. The quartic convergence of DMC4 is again verified. Its convergence is 
clearly distinct from lower order results and is nearly flat. In this case, we have no analytical 
results and all lines are just least square fits to the data. 

To demonstrate that DMC4 can be used to solve realistic physical problems, we use 
it to solve for ground state properties of bulk liquid helium described by the many-body 
Hamiltonian 

# = E-|^ + £^), (37) 



where % /m = 12. 12 A 2 K with potential V determined by Aziz et al [ 16| . Instead of the 
usual McMillan trial function, we use a trial function of the form, 

0(x) = JJexp{-lii(2) exp[-(r^ - c )/d ]}. (38) 

i<j 

With Co = 2.8A and do = 0.48A, this trial function gives a slightly better energy of 5.886(5) 
K/particle. Since the standard calculation details [17] are well known, we will just describe 



the results as summarized in Fig. 3. Again, the convergence of our 4th order algorithm is 
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clearly quartic. The extrapolated values are -7.114(2) K for our 4th order algorithm and 
-7.111(2) K for the second order algorithm DMC2a. Both are in agreement with Boronat 



and Casulleras's |18j second order DMC result of -7.121(10) K. Note that very large steps 
can be used with algorithm DMC4, roughly ten times as large as those of second order 
algorithms. 

In Fig. 4, we show the resulting radial density distribution g{r) from our 4th order cal- 
culation at e = 0.1 and e = 0.2. The distribution is virtually unchanged even at these 
large time steps and both are in excellent agreement with the experimental distribution of 
Svensson et al. [ |ITfl . 

In Fig. 5 we show DMC4's thermalization toward the exact ground state from the varia- 
tional trial wave function. Starting from the initial variational energy, we plot the population 
averaged energy as a function of iterated time for various time step sizes. This plot shows 
that each iteration of the algorithm at e = 0.2 is indistinguishable from multiple iterations 
at smaller time steps having the same time interval. Moreover, it demonstrates that the 
algorithm converges to the ground state inversely proportional to the step size used, up to 
e = 0.2. That is, only five iterations are needed at e = 0.2 to reach the exact ground state 
near t = 1.0 and 20 iterations at e = 0.05, etc.. Thus our 4th order algorithm can project 
out the ground state with ten times fewer updates than a second order algorithm. 

More important than the thermalization time is the observable correlation time. In a 
Monte Carlo calculation, it is highly desirable to have uncorrelated configurations for an 
accurate estimate of the statistical errors. The correlation coefficient for an observable O is 
defined by 

Co{At) = <o(t)o(t)>-<o(t)>* ■ m 

In Fig. 6, we show the ground state energy correlation function of liquid helium as computed 
by our 4th order algorithm. The correlation time is roughly At « 1.5, at which point the 
correlation coefficient dropped to zero. This plot shows that the correlation time depends 
only on the total time separation. Thus if the algorithm remains accurate at large time 
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steps, then fewer iterations are needed to produce uncorrected configurations. 

In implementing the fourth order Langevin algorithm, we used the standard fourth order 
Runge-Kutta algorithm to solve the trajectory equation (|T5|). When the step size is large, 
the 4th order error in the Runge-Kutta algorithm can greatly overshadow the intrinsic 4th 
order step size error of the Langevin and that of the DMC algorithm, causing both to fail 
prematurely. To guard against this from happening, we monitor the difference between 
the results of the fourth order Runge-Kutta and its embedded second order algorithm. If 
the square of this difference is larger than some tolerance, say 0.01 , we recalculate the 
trajectory twice at half the time step size. Even at the largest step size used, only a few 
percent of trajectories need to be recalculated, incurring only a small additional overhead. 
This additional effort greatly extended the flatness of the convergence curve as shown in 
Fig. 3. The total calculation time running on a single processor of an Origin 2000 is 51 
hours, about 5 times as long the second order algorithm. 

V. ALTERNATIVE FOURTH ORDER ALGORITHMS 



Our DMC4 algorithm (|23|) , which preserves the Fokker-Planck operator L intact, may 
not be the most efficient fourth order algorithm possible. Consider the limit in which the 
trial function approaches the exact ground state, 0(x) — > ipo(x). In this ideal case the local 
energy becomes an irrelevant constant, -El(x) — > Eq, and the algorithm is just 

which is the running of the 4th order Langevin algorithm twice, at half the time step. It 
seems plausible that one should be able to derive a 4th order DMC algorithm which reduces 
to a single run of the 4th order Langevin algorithm in the same limit. 

We are thus led to consider the general factorization, to fourth order, of a three-operator 
exponential e - e ( T + D + E i>) _ There are now 9 double commutators to be considered: 6 are the 
generalizations of the two operator case, 
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[T,[D,T]}, [D,[T,D]], [D,[E L ,D]], [E L ,[D,E L ]], [E L ,[T,E L }}, [T,[E L ,T}\, (41) 

and three new ones related by the Jacobi identity, 

[T, [D, E L ]] + [D, [E L , T]] + [E L , [T, D]] = 0. (42) 

Thus only two of the last three commutators are independent. Note also that for the present 
form of the operators, [El, [D,El\] = 0. We have examined all these commutators in the 
case of liquid helium to determine which one is doable and can be kept . To explore the many 
possible factorizations, we have devised a Mathematica program to combine the exponential 
of operators symbolically. With the help of this program, we have explored an extensive 
list of distinct 4th order algorithms. Since there are many operators in each such factor- 
ization, it is too cumbersome to write out the explicit exponential form. Moreover, since 
the factorization will always be left-right symmetric, there is no need to repeat operators on 
the left side. In the following, we will only indicate the exponential operators symbolically 
beginning with the central one and list only operators to the right. For example, algorithm 
DMC4 (|D will be denoted as 



T + D + E L ^ -E L +-L+-E L 



2 ~ 1 
3^ + 2 



-(1 i=)T + -D + —=f +-D + -(1 =)T 

2 K VS } 2 ^v/3 2 2 l v/3 ; 



+ \e l . (43) 



Each update of this algorithm requires the evaluation of, in decreasing order of computational 
complexity, 4D's, 2T's, IE L , \E L and 4T's. (The last E L from the last update can be used 
as first E L of the current update.) Since D is the most computationally intensive operator, 
followed by T, El, etc., we would like to minimize their appearance in that order. Below, 
we will describe two alternative algorithms that are computationally more economical than 
DMC4 in solving for the ground state of liquid helium. 

One possible 4th order algorithm is to retain the same double commutators [D, [T, D]] 
and [El, [T, El}} as in DMC4, but allow L = T + D to be broken up: 

T + D + E L « \f + \d + \e l + \t + l -D + l -T + \e l . (44) 

6 D 8 D 6 o 8 
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Here, T and E given by 



e 2 



T = T + -[D,[T,D]], (45) 

E L = £ L + ^[E L ,[T,£ L ]]. (46) 

This algorithm requires 4D's, but only IT, 1-El, 2E^s and 4T's. We will denote this algo- 
rithm as DMC4a. This algorithm is roughly 10% faster DMC4 and its quartic convergence is 
clearly demonstrated in Fig|| However, its convergence range is only about half of DMC4. 
The actual running time for this algorithm is 46 hours. 

To reduce the number of D operators, one must pay the price of retaining additional 
double commutators. We will refer to the following algorithm with only two D operators as 
DMC4b: 

T + D + E L « ^=a f + ^j=E L + ^=(1 - a )T + l -D + -^T + ~c E L + ±c T, (47) 



where a® = 1/y 1 + v^, Co = 1 — l/\/3 and 

«o 



-L«of = ^OoT+Ij 



(2 - V3)(\D, [T, £>]] + [£>, [B t , D]]) + (o - -|) [B t , [T, E L 



+ ^T. (48) 

The additional commutator 

[D, [E L: D]] = -GiViiGjVjELfr)], (49) 

is a calculable function involving higher derivative of Gi and El- In this algorithm, we 
have placed all the double commutators at the center so that they are evaluated only once 
per update. This is done by splitting -^aoT ~~ * ^m a oT + ...j^a T in Eq. (f£8|) , meaning 
that we first do half of required Gaussian random walk, evaluate all double commutators, 
then complete the remaining half of the Gaussian walk including the effect of [D, [T, D] 
as it is done in the Langevin algorithm. The ubiquitous irrational coefficients are roots of 
quadratic equations which force unwanted double commutators to vanish. One can check 
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by inspection that as 0(x) — ► ipo(x) and -Ex(x) — > E , ([47]) reduces to just the 4th order 



Langevin algorithm ( |26| ) 



The results of this algorithm for liquid helium are shown in Fig.[| as asterisks. The 
ground state energy is correctly obtained, but because these higher derivatives are rapidly 
varying functions, we have not been able to stabilize the population of weights beyond 
e ~ 0.05. While this algorithm may not work as well as DMC4 and DMC4a for liquid 
helium, its economy of requiring only two trajectories per update may be of utility in other 
applications. The calculation time for data points shown is 36 hours. 

VI. CONCLUSIONS 

In this work, we have derived a number of distinct fourth order DMC algorithms by 
factorizing the imaginary time Schrodinger evolution operator to fourth order. This is a 
notable advance in algorithm development, made possible only by the recent progress in 
understanding positive coefficient operator factorization. Our work illustrates a global view 
of algorithms as products of factorized operators. Such a perspective gives order and insight 
into the working of the Diffusion Monte Carlo algorithm. It would have been very difficult 
to derive such a high order algorithm without such a conceptual structure. We have further 
demonstrated the practicality of these algorithms by using them to solve for the ground 
state of liquid helium. The quartic convergences of DMC4 and DMC4a have been verified 
and both yielded ground state energy and radial density distribution in excellent agreements 
with experiment. Despite the fact that these algorithms are rather complicated to program 
requiring higher order derivatives, they allow very large step sizes to be used, virtually elim- 
inate the time step size error and greatly reduce statistical correlations between successive 
updated configurations. 

Since this is only the first demonstration of quartic algorithms, there is room for further 
improvements. For example, we have shown how the factorization of a three-operator ex- 
ponential can lead to a number of distinct four order DMC algorithms. A more systematic 
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categorization of various fourth order factorizations would help in obtaining the most efficient 
algorithm. Secondly, the retainment of some double commutators is necessary, however, it 
has not been studied in detail where they should be placed so as to minimize the 4th order 
error coefficient or computational effort. It is observed that the step size convergence curve 
is flatter when the trajectory equation is solved more exactly. In this work, we have only 
used the 4th order Runge-Kutta algorithm in solving for the deterministic trajectory. Future 
study may explore the effects of using alternative numerical methods, such as symplectic 
algorithms ||20|| , for solving the trajectory equation. 
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FIGURES 
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FIG. 1. The convergence of the ground state energy of a 3-D harmonic oscillator as a function 
of the time step size e. The diamonds are first order DMC1 simulation results. The filled and open 
triangles are second order DMC2a and DMC2b results. The asterisks indicate results of a linear 
algorithm with rejection. The filled circles are DMC4 results, Eq.(23). The errorbars are smaller 
than the symbol size. The various lines are the corresponding exact analytical results except in 
the case of the rejection algorithm. For the latter case the line is just a 6th order polynomial fit to 
the simulation data. 
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FIG. 2. The convergence of the ground state energy of a 3-D Morse oscillator as a function of 
the time step size e. The diamonds are first order DMC1 simulation results. The filled and open 
triangles are second order DMC2a and DMC2b results. The filled circles are the 4th order results 
of DMC4. The various lines are least square fits to the simulation data. The first order results are 
fitted with a parabola, the second order results by a cubic polynomial, and the fourth order results 
by just a constant plus a fourth order term in e. 
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FIG. 3. The time step convergence of the ground state energy per particle for bulk liquid helium 
in a 128 particle simulation. The solid circles are the result of our 4th order algorithm DMC4. The 
open circles and asterisk are for algorithm DMC4a and DMC4b respectively. For comparison, we 
also show as triangles, second order results from algorithm DMC2a. The lines are least square fits 
to the data. 
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FIG. 4. The radial density distribution of bulk liquid helium. The circles are DMC4 results at 
e = 0.1 and the crosses are DMC4 results at e = 0.2. The solid line is the experimentally extracted 
g{r) of Svensson et al. (Til at 1 K. 
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FIG. 5. The relaxation of liquid helium's ground state energy toward its exact value as sim- 
ulated by DMC4 at three time step sizes. The large circles are at e = 0.2 and the small circles 
are at e = 0.05. They are connected by straight line segments to guide the eye. The dotted line 
corresponds to results at e = 0.01. 
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FIG. 6. The correlation coefficient function, Eq.(^), for the ground state energy of liquid 
helium as computed by DMC4 at three time step size of e = 0.2 (large circles), e = 0.05 (smaller 
circles) and e = 0.01 (dotted line). The connecting line segments are for guiding the eye only. 
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